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DESCRIPTION 

METHOD OF COMPUTING FIR FILTER COEFFICIENT AND PROGRAM 

FOR COMPUTING SAME 

CROSS REFERENCE TO RELATED APPLICATION 

This application is based upon and claims the benefit of priority from the prior 
U.S. provisional Patent Application No. 60/418,313, filed on October 15, 2002, the 
entire contents of which are incorporated herein by reference, 

TECHNICAL FIELD 

The invention relates to a computational method for impulse response 
coefficients of a universal maximally flat FIR filter, and to a computational program for 
the same. 

BACKGROUND ART 

More than thirty years have passed since the publication of the seminal paper of 
p. Herrmann on the design of maximally flat FIR digital filters [1]. Of the body of the 
related research that has been published since Herrmann's article, the most significant 
contribution to the theory of lowpass maximally flat FIR filters has been Baher*s results 
on the design of filters with simultaneous conditions on. amplitude and group-delay 
response [2]. In recent years, it has been recognized that Baher's filters form a unifying 
class of maximally flat filters. This class encompasses various types of seemingly 
distinct systems [3]. Calling this class of filters the universal maximally flat digital 
filters, the authors simplified Baher's formula for the transfer function and showed that 
lowpass filters of even or odd lengths, with linear or nonlinear phase response 
characteristics, as well as fractional delay and half-band filters may be readily obtained 
by setting the design parameters in an. appropriate manner. 

Notwithstanding universal maximally flat digital filters have an explicit formula 
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for their transfer functions and, thus, may be readily designed using a computer algebra 
system, the simplified formulas provided in [3] are still unwieldy. A direct use of the 
formulas involves computation of binomial coefficients with possibly non-integer 
indices that appear as terms of the summand in a three-fold sum. Binomial coefficients 
5 with non-integer indices cannot be computed by the simple method of Pascal's triangle. 
Another common problem with computations involving the binomial coefficients is due 
to their tendency to introduce very large integers to the intermediate steps of the 
computation. With the flexible computer power today, this may not be a serious 
impediment in many situations. However, there can be cases where limitations on 
10 hardware and software resources dictate the use of efficient means in computation of 
filter coefficients. It is desirable that such methods be free of binomial coefficients and 
have a simple iterative structure. Moreover, it is known that if the design parameter 
associated with the group delay is limited to rational numbers, as is the case with many 
realistic settings, the values of the impulse response coefficients are rational numbers. In 
15 such cases, it is important to ensure that the answer to the numerical problem in 
computation of an impulse response coefficient is exactly a fraction, say 1/3, not a 
floating number that may be shown as "0.333333335". Therefore, it is desirable for a 

computation algorithm to maintain this property of the impulse response coefficients 

and yield rational values for rational design parameters. 
20 The following are prior art references related to the present invention. 

[1] O. Herrmann, "On the approximation problem in nonrecursive digital filter design," 

IEEE Trans. Circuit Theory, vol.CT-18, no.3, pp.411-413. May 1971. 

[2] H. Baher, "FIR digital filters with simultaneous conditions on amplitude and delay," 

Electron. Lett., vol. 18, pp296-297, 1982. 
25 [3] S. Samadi, A. Nishihara and H. Iwakura, "Universal maximally flat lowpass FIR 

systems," IEEE Trans. Signal Processing, vol.48, pp. 1956-1964, 2000. 

[4] R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics. 

Addison-Wesley, 1989.' 



DISCLOSTTRE OF INVENTION 

An object of invention is to compute impulse response coefficients of the 
universal maximally flat FIR filter efficiently in a short period of time. 

According to one aspect of the present invention, the computer executes a first 
operation by a first recurrence formula, receiving a filter order (positive integer) of a 
universal maximally flat FIR filter, the number of zeros at z=-l (integer equal to or more 
than zero), and a parameter for a group delay at z=l (rational number). The recurrence 
formula includes parameters the filter order, the number of zeros, and the group delay, 
and provides coefficients in Bernstein form representation of a transfer function of a 
universal maximally flat FIR filter. 

The first recurrence formula, for example, is expressed as 

bj' = (-l){(2d) bj.i' + ( j - 1 ) bj.2'} / (N - j + 1 ) where 1 ^j^N with bo' = 1 
andb-i' = 0, 

in which the filter order is N, the parameter for the group delay is d, the coefficients in 
Bernstein form representation of a transfer function of the universal maximally flat FIR 
filter are bj'. The resultant of the first operation is expressed as 
B'=|^l^br,...,bN.K',0,...,0} where the number of zeros is K. 

Using a resultant of the first operation as an initial value, the computer executes 
the second operation by the second recurrence formula composed of additions, 
subtractions, and divisions by 2 to extract impulse response coefficients of the 
universal maximally flat FIR filter. 

The second recurrence formula, for example, is expressed as 
hi^P^ = ( 1 + E ) hi^P-^^ / 2 + (1 - E) hi.i^P-^^ / 2 where 1 ^p^N, O^i^p with ho^°^ 
= B' and h./P>= {0,...,0}, 
in which a sequence for computing impulse response coefficients of the universal 
maximally flat FIR filter is expressed as hi^P^=(hij^P^)=(hi,o^P\hia^P\...). and an arbitrary 
sequence Ai is expressed as E^ = E ( E-^'^Ai )^ E^Aj = EAi = Ai+i^ E°Ai = Ai in which a 
forward shift operator satisfying the expression is E. The impulse response 
coefficients extracted from the resultant of the second operation are expressed as 
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hi=hi/^^ where O^i^N. 

RRTEF DESCRIPTION OF DRAWINGS 

Fig. 1 illustrates steps of introducing a computational algorithm for impulse 
5 response coefficients according to the present invention. 

Fig. 2 illustrates steps of proving the proposition 1. 
Fig. 3 illustrates steps of proving the corollary 1. 

Fig. 4 illustrates an overviev^ of the computational algorithm for the impulse 
response coefficients according to the present invention. 
10 Fig. 5 is a flowchart showing the computational algorithm for the impulse 

response coefficients. 

Fig. 6 is a program list of the computational algorithm for the impulse response, 
coefficients executed by computer. 

15 BEST MODE FOR CARRYING OUT THE INVENTION 

The embodiments of the present invention will be described with reference to 
the drawings in the following. 

First, a description will be made on mathematical underpinning of the 
computational algorithm for the impulse response coefficients according to the present 
20 invention. Fig. 1 illustrates steps of introducing the computational algorithm for the 
impulse response coefficients. 

It was shown in [3] that the transfer function of the universal maximally flat 
digital filters, a family of filters identical to those proposed by Baher [2], is given by 

HN.K.dW= H b,(l^)'(l:^)--, (1) 

0<><N— K 

25 where the coefficients bj are characterized by the three-terni recurrence 

j bj + 2d bj«i - 0 N - 2) b^_2 = 0, j > 1 , (2) 



with the initial values 



bo = 1 and b_i = 0. (3) 

By converting (1) to the power form representation 

HN,K.dW= Hicz-^ (4) 

0<k<N 

the impulse response coefficients hk become [3] 

^ 0<j<N-K 0<p<lcO<i<j 

Obviously, (5) is a computationally expensive formula. To ameliorate this 
situation, devised is an algorithm that exploits the combination of (1) and (2) to compute 
hkfork = 0,...,N. 

First, (1) is written in the traditional Bernstein form by introducing the sequence 
10 bj' defined by 

I 0<j<N, 
hl^^'l (?) (6) 



0, otherwise. 



Then obtained is 



0<j<N-K \ ' / 



0<j<N- 

This is the Bernstein form of the transfer function in the traditional sense. 
15 Although (7) has an additional binomial coefficient compared to (1), we will see later 
that the introduction of bj' results in considerable simplification. 

A recurrence for bj' is obtained by substituting (6) into (2) and dividing both 
sides by j (^) ^0. 

Thus, 

b/ + 2d^b|_,-^{i-N-2)^b;_, = 0, 1<i<N. (8) 
20 ^^'^ 

The above recurrence can be simplified by invoking the definition of the 

binomial coefficients [4] and a few simple algebraic steps. It follows that 
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with the initial values 

bA = l and hL^^O. (10) 

Thus, for given values of N and d the values of bj' can be computed recursively by (9). 
5 The next step is to convert (7) to the power form (4) and obtain the impulse 

response coefficients. A direct expansion of the summand in (7) using the binomial 
theorem will introduce additional sums and new binomial coefficients. Hence, this is not 
a desirable action to take. The other alternative is to use a conversion identity well 
known in the mathematics literature and in the field of computer-aided geometric design. 
10 The identity is based on taking successive differences of the Bernstein coefficients. We 
state this identity as a proposition and provide a proof for it. 
(Proposition 1) 

Allow the Bernstein form of the polynomial p(x) = Z o^i^N Pi x' to be given by 

p(x)= wMxMi-x)-'-^ (11) 

0<i<N N ^ / 

15 Then 

.^bo, i = 0 N, (12) 



where the operator A is the forward difference operator. 

The forward difference operator is defined by the geneiral recurrence 

A^K = A(A^-^bO, (13) 

20 together with 

A^bi= Abi = bi+i -W, A^W = bi. (14) 

For instance, A bo = bi - bo and A.^bo = b2 - 2bi + be In this notation, the 
operator acts on the indices, written as subscripts, of the members of a sequence. To 
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prove the proposition, we also need the forward shift operator E. The effect of E on the 
index i of a sequence like bi is specified by 

E'bi= E(E-^bO. (15) 

together with 

bi = Ebi = bi+i . E° bt = bi. (16) 

5 

E is related to A by 

A = E-1. (17) 

(Proof of proposition 1) 

Fig. 2 illustrates steps of proving the proposition 1, Starting with the Bernstein 
10 form of p(x) and expand the term (1 - x)^"* using the binomial theorem to obtain 

0<i<NO<j<N-i \ / \ ^ / 

By the symmetry of the binomial coefficients and the trinomial revision 
property [4], obtained is 

a)rr^(^)(.-Ti,)- 

15 On the other hand 

bt = E^bo. (20) 

Thus, it can be written as 

0<v<N 0<j<N-t \ ' / \ '/ 

The ranges for the indices of the two-fold sum above may be interchanged as 
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Z Z = E L (22) 

0<i<N 0<j<N~i 0<i<N 0<i<j 

with no effect on the value of the sum. Now, N - j is substituted for j in the sum, and 
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obtained is 



0<N-j<NO<i<j V V \ V 



(23) 



The right side above simplifies to 



(25) 



(26) 



0<j<N 0<i<j 

5 Note that 

H (-1 y-' (. \) r bo = ( t - 1 )^ bo. 

0<l<j V V 

It then follows that 

pw= 21 ^^f^)(t-iybo 

0<j<N ^ ^ / 

that is the power form representation of p(x). From the uniqueness of the power from 
10 representation of a polynomial and the relation between the operators E and A it 
follows that Pi = (^) A 'bo. (the end of the proof) 

Now presented is a main and last result in this preparatory phase before 
proceeding to the derivation of the algorithm. 
(Corollary 1) 

15 The transfer function HN.K,d(z) is given by 

HN.K,d(^) = + )^ K (27) 

(Proof of corollary 1) 

Fig. 3 illustrates steps of proving the corollary 1. 
Let 

= l-2x. (28) 

20 

Then, it can be written that 



HN.K.dW= H W(^)xi(l-x)^^-'. (29) 

0<j<N-K \ ^ / 

By the above proposition obtained is 

HN.K.d(x)= H (^)A'bix^ (30) 

0<j<N \ ^ J 

By substituting for x in terms of z'^ in the above relation and expanding the 
summand using the binomial theorem, obtained is 

By the trinomial revision identity obtained is 

The sum indices can also be modified to write 

H.,K..(.) = L (^) G^f i^^-K. (33) 

0<i<N \ / i<j \^ / 

Using p = j - i as an index for the second sum, obtained is 

0<i<N \ / 0<p<N-i \ ^ / 

This simplifies to 

Hn.k..(z)= 21 (^)(::^)Mi+|)^-X. (35) 

0<i<N \ / 

Invoking the binomial theorem further makes it possible to write 

HN,K,dW = (l+f-f^-M^bi.' (36) 

By using the relation between A and E, (27) follows, (the end of the proof) 

Next a description will be made on details of the computational algorithm for 

the impulse response coefficients according to the present invention. 

The algorithm consists of two major stages. In the first stage, by running the 

recurrence (9) from j = 1 up to j = N - K, the following sequence is obtained. 
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B' = (1,b{...,bj^_K,0....,0). (37) 

The entries bj' of B' for N - K < j ^ N are set to zero. The next step is to 
expand the transfer function using (27). This can be done by noting that 

5 for arbitrary integer p > 0. Evaluation of both sides above in the power form yields 

^h(p),-< = (1+1 + l^z-i) ^U^-''z-\ (39) 

where hi^^^ is a sequence of coefficients defined by 

hl^^ = (hi^?)=(l4:o^Hg^••o> (40) 

and operator E effects a forward shift on the index j. Thus, it is obtained that 

Y_ hl-^z-^ = ^ h1^-^5 z-^ + 1^ Hi--^^ z-^-^) . (41) 

10 1 i 

It follows that 

HjPl = I+IhI^-^J- 4- 1^h1IT^\ (42) 

The ranges for i and p are specified by 

l<p<N, 0<i<p. (43) 

15 The initial values are given by 

h!^^^B', H!r,^ = (0.0....). (44) 

This completely specifies the recurrence we use in the second stage of our 
algorithm. After completing the Nth step, at the point when p has reached N, the value 
of the i-th impulse response coefficient is given by 

Hi = first term of sequence h['^\ i = 0, . . . , N . (45) 
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Fig. 4 illustrates an overview of the computational algorithm for the impulse 
response coefficients according to the present invention. In Fig. 4, each node represents 
a sequence hi^''^ The nodes marked by O's represent zero sequences. The recurrence is 
illustrated in Fig.4 for N = 3 and an unspecified value of K. It can be visually verified 
5 that the recurrence has a triangular structure. 

Fig. 5 is a flowchart showing the computational algorithm for the impulse 
response coefficients according to the present invention. 

Receiving a filter order N (positive integer) of the universal maximally flat FIR 
filter, a number of zeros K (integer equal to or more than zero) at z=-l, and a parameter 
10 d (rational number) for a group delay at z=l, the computer computes impulse response 
coefficients according to the following steps. 

In step SIO, all of N+1 elements B'[j] (O^j^N) of a single-dimension array B' 
are initialized to zero. In step S20, all of elements r[p,ij] (0 ^ p ^ N, 0 ^ i ^ N, 0 ^ j 
^N) of a three-dimension array r are initialized to zero. In step S30, an element B'[0] 
15 of the single-dimension array B' is set to 1. 

In step S40, every element of the single-dimension array B' is decided by 
changing in sequence an index j from 1 to N-K in a recurrence formula B'|j] = (-l)X 
{(2d)B'[j-l] + 0-l)B'D-2]} / (N - j + 1). That is, the operation by the recurrence formula 
(9) is executed. 

20 In step S50, elements r[0,0,j] (0 ^N-K) of the three-dimension array r are set 

to elements B'|j] (O^j ^N-K) of the single-dimension array B'. 

In step S60, every element of the three-dimension array r is decided by 
sequentially changing, in the order of indexes j, i, p, an index j from zero to N-p, and an 
index i from zero to p, an index p from 1 to N in a recurrence formula r[p,i,j] = 
25 ( r[p-l,i-l,j] - r[p-l,i-l,j+l] ) / 2 + ( r[p-l,i,j] + r[p-l,i,j+l] ) / 2. That is, the operation 
by the recurrence formula (42) is executed. 

In step S70, N+1 elements h[i] (O^i^N) of a single-dimension array h are set to 
elements r[N,i,0] (O^i^N) of the three-dimension array r. 

In step S80, the single-dimension array h is output as impulse response 
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coefficients of the universal maximally flat FIR filter. This completes the 
computational processing of the impulse response coefficients by the computer. 

Fig. 6 illustrates a program list of the computational algorithm for the impulse 
response coefficients. by computer. A pseudo-code of a procedure that returns the value 
5 of impulse response coefficients as an array is given in Fig. 6, 

The following example helps elucidate the details. Let N = 3, K = 1 and d = 
(_l)/4. The actual values of the sequences hi^^^ i = 0,...,3, corresponding to the 3-D 
array r in the procedure of Fig. 6, are as follows: 

h(p) hi-^ h^^ K^' 

p=0: (1.1/6,-11/24,0) 

p^l. (7/12,-7/48.-11/48) (5/12,5/16,-11/48) (4^) 

p=2: (7/32,-3/16) (35/48,1/12) (5/96.13/48) 

p=3: (1/64) (39/64) (31/64) (-7/64) 

10 Only non-zero values of the sequences hi^^^ , i = 1,...,3, are presented. Note that 

all computations have been done in the field of rational numbers in an exact manner. 
Hence, round-off errors are avoided. 

As described above, an alternative representation for the transfer function of 
universal maximally flat FIR filters enables the development of a very simple algorithm 

15 for evaluation of the impulse response coefficients of the filters. The alternative form of 
the transfer function is based on the shift operator E and is derived from a technique 
referred to as finite calculus in [4]. 

The algorithm consists of two stages. In the first stage a two-term recurrence 
- relation is used to obtain a 1-D sequence of numbers. If the value of the group delay 

20 parameter d is a rational number, then the outputs of the recurrence are rational. The 
outputs are then used in a simple N-step algorithm of triangular structure. In each step 
of this algorithm, the values of the sequences obtained in the preceding step undergo 
simple additions, subtractions and divisions by 2 to yield sequences that are shorter in 
size. Upon the completion of the Nth step, N sequences are obtained. The first entries of 
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these sequences are the values of the impulse response coefficients. As long as the 
computations in the second stage are done in the field of rational numbers, the final 
values remain rational. This provides an accurate and fast means for exact computation 
of the impulse response coefficients. 

Here, a worst-case analysis of the computational complexity of the present 
algorithm will be described. 

The computational complexity may be analyzed by evaluating the complexity of 
the simple for-loop marked by (A) and the nested for-loop designated by (B) in Fig.6. 
For the number of additions and subtractions C±, it is obtained that 

C±=^4+2;^J2 Z^. (47) 

1<j<N-K 1<P<N 0<-L<p 0<j<N-p 

For a closed-form expression, note that 

C± = 4(1S|-K:)H-3 ^ Y_ (N-p + 1), . (48) 

l<P<NO<i<p-hl 

that can be evaluated as 

C±=4(N-K:) + 3 ^ (pH-2)(N-p). (49) 

0<p<N 

Thus, obtained is 

C± = 4(N-K)-h3{|N+N^ + ^N^). (50) 

For the number of multiplications C* obtained is 

C*= X 2 = 2(1SI-K). (51) 

l<j<N-K 

The number of divisions C/ is given by 

c/= r 1+ r z i: 2. (52) 

l<j<N-K l<p<N 0<i<pO<j<N-p 

that becomes 

C/ = (N - K) -h 2(| N -h + ^ N^). (53) 
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In comparison, for the direct evaluation using (5), the number of additions and 



subtractions 0+"^'"" is 



cDirect> ^ ^ 

0<j<N-K 0<p<kO<i<j 



The right side above is of the order OCN"*). The number of multiplications C* '^^"^ , on the 
5 other hand, is 

^Direct > ^ L 21 ^ ^- (^^^ 

0<j<N-KO<p<kO<i<i 

The right side above is of the order 0(N^). For the number of divisions C/^''^''^ in a direct 
evaluation, it can be written in 

0<j<N-KO<p<kO<i<j 

10 The right side above, again, is of order 0(N^). Thus, in addition to the simplicity of 
most of the division operations in the present algorithm (that are divisions by 2), a 
significant reduction in the overall complexity is accommodated in terms of sheer 
numbers. 

The invention is not limited to the above embodiments and various 
15 modifications may be made without departing from the spirit and scope of the invention. 
Any improvement may be made in part or all of the components. 

INDUSTRIAL APPLICABILITY 

As described above, the impulse response coefficients of the universal 
20 maximally flat FIR filters may be computed efficiently using a two-stage algorithm. For 
a filter of order N with K zeros at z = -1, the first stage of the algorithm uses a two-term 
recurrence to generate a sequence of length N - K -f 1. This sequence undergoes an 
N-step iterative process that involves additions, subtractions and divisions by 2. The 
algorithm does not involve any form of binomial coefficients. Moreover, the algorithm 
25 requires less number of arithmetical operations compared to a direct evaluation using 
the closed-form formula. The algorithm is mostly suitable for real-time generation of the 
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coefficients on a DSP chip because it may be easily implemented in a computational 
environment that has restrictions on the available hardware or software resources. 
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CLAIMS 

WHAT IS CLAIMED IS: 

1. A method of computing FIR filter coefficients, comprising the steps of.. 

inputting a filter order of a universal maximally flat FIR filter, a number of 
zeros at z=-l, and a parameter for a group delay at z=l, the filter order being a positive 
integer, the number of zeros being an integer equal to or more than zero, the parameter 
being a rational number; 

executing a first operation by a first recurrence formula which, includes 
parameters for the filter order, the number of zeros, and the group delay, and provides 
coefficients in Bernstein form representation of a transfer function of the universal 
maximally flat FIR filter; 

executing a second operation by a second recurrence formula composed of 
additions, subtractions, and divisions by 2, by using a resultant of the first operation as 
an initial value; and 

extracting impulse response coefficients of the universal maximally flat FIR 
filter from a resultant of the second operation. 

2. The method according to claim 1, wherein: 

the first recurrence formula is expressed as 

bj' = (-l){(2d) bj-i' + ( j - 1 ) bj-2'} / (N - j + 1 ) where 1 ^j^N with bo' = 1 
and b-i' = 0, 

wherein the filter order is N, the parameter for the group delay is d, coefficients in 
Bernstein form representation of a transfer function of the universal maximally flat FIR 
filter are bj'; 

the resultant of the first operation is expressed as B'={l,bi',...,bN-K',0,...,0}, 
wherein the number of zeros is K; 

. the second recurrence formula is expressed as 

hi^P^ = ( 1+ E ) hi^P-^^ / 2 + (1 - E) hi./P-^^ / 2 where 1 ^p^N, O^i^p with ho^°^ 
= B' andh./''^= {0,...,0}, 
wherein a sequence for computing impulse response coefficients of the universal 
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maximally flat FIR filter is expressed as hi^''HK^''^)Hho^^\Ki^''\ ■ ■ X and an arbitrary 
sequence Ai is expressed as = E ( E^'^Ai ), E'A; = EAi = Ai+i.E°Ai = Aj in which a 
forward shift operator satisfying the expression is E; and 

the impulse response coefficients extracted from the resultant of the second 
operation are expressed as hi = hj.o^^^ where O^i^N 

3 A program for computing FIR filter coefficients, the program causing a computer to 
execute the steps of: 

determining every element of a single-dimension array B' using a filter order N 
being a positive integer of a universal maximally flat FIR filter, a number of zeros K at 
z=-l, K being an integer equal to or more than zero, and a parameter d for a group delay 
at z=l, d being a rational number, all of which are provided by inputs, by changing in 
sequence an index j from 1 to N-K in a recurrence formula B'[j] = (-1) X {(2d)B'|j-l] + 
(j-l)B'[j-2]} / (N - j + 1), the single-dimension array having N+1 elements B'jj] where 0 
^N, in which an element B'[0] thereof is initialized to 1 and all the elements thereof 
except the element B'[0] are initialized to zero; 

determining every element of a three-dimension array r by sequentially 
changing, in the order of indexes j, i, p, an index j from 0 to N-p, and an index i from 0 
to p, an index p from 1 to N in a recurrence formula r[p,i,j] = ( r[p-l,i-l,j] - 
r[p-l,i-l,j+l] ) / 2 + ( r[p-l,i,j] + r[p-l,i,j+l] ) / 2, the three-dimension array r having 
elements r[p,i,j] where O^p^N, O^i^N, O^j^N, in which elements r[0,0,j] thereof 
where O^j ^N-K are initialized to elements of the single-dimension array B'[j] where 0 
^j^N-K, and all the elements thereof except the elements r[0,0,j] are initialized to 
zero; and 

extracting elements r[N,i,0] of the three-dimension array .r where O^i^N as 
the impulse response coefficients of the universal maximally flat FIR filter 
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ABSTRACT 

The computer executes a first operation by a first recurrence formula, receiving 
a filter order (positive integer) of a universal maximally flat FIR filter, the number of 
5 zeros at z=-l (integer equal to or more than zero), and a parameter for a group delay at 
z=l (rational number). The first recurrence formula includes parameters for the filter 
order, the number of zeros, and the group delay, and provides coefficients in Bernstein 
form representation of a transfer function of a universal maximally flat FIR filter. The 
computer then executes a second operation composed of additions, subtractions, and 
10 division by 2 by a second recurrence formula by using a resultant of the first operation 
as an initial value to extract impulse response coefficients of the universal maximally 
flat FIR filter from a resultant of the second operation. 
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START 



set B' [j]=0 for j=0 to N 



set r [p, i , j]=0 for j. i . p=0 to N 



set B' [0]=1 



SIO 



S20 



S30 



run recurrence 
8* [j]=(-l) {(2d)B' [j-1] + (j-l)B' [j-2]l/(N-j+1) 
for j=l to N-K 



set r[0,0, j]=B' [j] for j=0 to N-K 



run recurrence 
r[p. i. j]=(r[p-1. j]-r[p-l, i-1. j+1])/2 
+ (r[p-1. i. j]+r[p-l, i. j+1])/2 
for j=0 to N-p. i=0 to p, p=1 to N 



set h[i]=r[N. i.O] for i=0 to N 



output h[i] for i=0 to N 



S40 



S50 



S60 



S70 



S80 



END 



Fig. 5 
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